R as a GIS
2019/5/7
1 はじめに
RをGISとして使う． ポリゴンオブジェクト上に位置するラインオブジェクトの延長を算出する． それぞれ以下のようにする．
- ポリゴン：土砂災害警戒区域
- ライン：道路
土砂災害警戒区域上に位置する道路を特定し，その延長を計算することで防災政策を立案するうえでの一つの判断材料を与える．
2 前準備
使用するライブラリの読み込み，国土数値情報からのダウンロード，openstreetmapからのダウンロードを行う．
2.1 ライブラリの読み込み
library(sf)
library(tidyverse)
library(tmap)
library(mapview)2.2 国土数値情報から
石川県のものについて，以下を読み込む．
- 行政区域：N03~
- 土砂災害警戒区域：A33-17~
コードは以下のとおりである．
# 日本語を扱うのでCP932でのエンコードを忘れない
isikawa <- st_read("N03-170101_17_GML/N03-17_17_170101.shp", options = "ENCODING=CP932", stringsAsFactors=F)
# emergency <- st_read("N10-15_17_GML/N10-15_17.shp", options = "ENCODING=CP932", stringsAsFactors=F) 不使用
sediment <- st_read("A33-17_17_GML/A33-17_17_GML/A33-17_17Polygon.shp", options = "ENCODING=CP932", stringsAsFactors=F)N03_004が市町村名に相当する．
ここを七尾市に指定すればよさそうである．
nanao <- isikawa %>%
filter(N03_004 == "七尾市")行政区域を七尾市のみに絞ったものができた． このポリゴン上で土砂災害警戒区域を絞り込む．
#emergency_nanao <- st_intersection(nanao, emergency)
sediment_nanao <- st_intersection(nanao, sediment)そして投影座標系をEPSG2449とする．
nanao <- st_transform(nanao, 2449)
#emergency_nanao <- st_transform(emergency_nanao, 2449)
sediment_nanao <- st_transform(sediment_nanao, 2449)mapview()を使ってインタラクティブに確認する．
mapview(nanao) +
mapview(sediment_nanao)2.3 openstreetmapから
openstreetmapのapiを使用して七尾市の道路情報を取得する． このapiを使えば指定した矩形範囲内の道路情報を得ることができる． ただし対象とする範囲が広すぎる場合は取得まで時間がかかりすぎることがある． 今回は国道や市道などの種別ごとに取得することで，一度に得るオブジェクトの数を減らすことで対処する． 国道や市道などの属性は下のような名称となる．
- morotway
- trunk
- primary
- secondary
- tertialy
- unclassified
- residential
矩形範囲は緯度経度での指定となる．
七尾市を囲む範囲は(36.964677, 136.772596,37.198499, 137.070257)のようになる．
openstreetmapのapiについてはpythonでのラッパーがあるのでそれを使用する．
Rからpythonを呼び出すにはreticulateパッケージをつかう．
pythonを使う下準備は次の通りとなる．
reticulate::use_python(Sys.which(names = "python"))
opass <- reticulate::import(module = "overpass")
api = opass$API() # Rからは$で潜っていく．たくさん作るので，関数にしておく
geojson_as_sf <- function(way_attr){
# way_attr は文字列
bound_box <- "36.964677, 136.772596,37.198499, 137.070257"
query_opass <- str_c('way["highway"="',way_attr,'"](',bound_box,');(._;>;)')
attr_api <- api$get(query_opass, verbosity = 'geom')
way_sf_NodeWay <- st_as_sf(st_read(attr_api))
way_sf <- way_sf_NodeWay %>%
filter(highway == way_attr)
way_sf <- st_transform(way_sf, 2449) # 4326から2449にする
}道路属性ごとに異なるsfオブジェクトを生成する．
道路属性を文字列として格納しgeojson -> sf する．
# やたらと出力が長くなる． message = F だとhtml出力が止まる．
attr_str <- c("motorway", "trunk", "primary", "secondary", "tertiary", "unclassified", "residential")
for (i in attr_str) {
now_sf <- str_c(i, "_sf")
assign(now_sf, geojson_as_sf(i))
}バウンディングボックスは矩形であり，七尾市の外側の道路も取得している． これを七尾市上のみに絞り込む．
for (i in attr_str) {
now_sf <- str_c(i, "_sf")
nanao_sf <- str_c(i, "_nanao") # 道路のsfについて 七尾市だけのものをこれに格納する
assign(nanao_sf, st_intersection(nanao, get(now_sf))) # nanao_sf <- st_intersection(nanao, atr_sf) とするのをまあ
}最終的には全道路属性が1つのデータフレームにまとまっているsfオブジェクトを生成する．
現時点ではまだ，道路属性ごとに異なるsfオブジェクトである．
これを以下の手順でひとまとめにする．
- 道路属性個数分だけある
sfオブジェクトをデータフレームとsfcオブジェクトに分解． - 1.で，道路属性個数分のデータフレームができたので，それをひとまとめにする．
- 1.で，道路属性個数分の
sfcオブジェクトができたので，それをひとまとめにする． sfをつくる．2.でできた1つのデータフレームと，3.でできた1つのsfcをくっつける
分解作業，すぐおわる．
for(i in attr_str){
now_nanao <- str_c(i, "_nanao")
# データフレームだけを now_df に格納
now_df <- str_c(i, "_df")
assign(now_df, st_set_geometry(get(now_nanao), NULL))
# sfc だけを now_sfc に格納
now_sfc <- str_c(i, "_sfc")
assign(now_sfc, st_geometry(get(now_nanao)))
}bind_rowにより，道路属性ごとのデータフレームを下に付け足す．
# 格納用のデータフレームを作って，それにバインドしていく
Road_attr_df <- data.frame()
for(i in attr_str){
now_df <- get(str_c(i, "_df")) # ここでget()すれば1行でデータフレームになる
Road_attr_df <- bind_rows(Road_attr_df, now_df)
}sfcも全道路属性をひとまとめにする．
リストに格納する．
# 空のリストに放り込む
Road_sfc <- list()
for (i in attr_str) {
now_sfc <- get(str_c(i, "_sfc")) # ここでget()すれば1行で sfc
Road_sfc <- c(Road_sfc, i = now_sfc)
}sfオブジェクトをつくる．
Road_sf <- st_sf(Road_attr_df, geometry = Road_sfc)
Road_sf <- st_set_crs(Road_sf, 2449) # 本当は上のやつでいっしょにしたい道路属性ごとに色分けして図化する．
mapview(nanao) +
mapview(Road_sf,
zcol = "highway")